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Abstract 

Background: Habitat fragmentation has accelerated within the last century, but may have been ongoing over 
longer time scales. We analyzed the timing and genetic consequences of fragmentation in two isolated lake-dwelling 
brown trout populations. They are from the same river system (the Gudena River, Denmark) and have been isolated 
from downstream anadromous trout by dams established ca. 600-800 years ago. For reference, we included ten other 
anadromous populations and two hatchery strains. Based on analysis of 44 microsatellite loci we investigated if 
the lake populations have been naturally genetically differentiated from anadromous trout for thousands of years, 
or have diverged recently due to the establishment of dams. 

Results: Divergence time estimates were based on 1) Approximate Bayesian Computation and 2) a coalescent-based 
isolation-with-gene-flow model. Both methods suggested divergence times ca. 600-800 years bp, providing strong 
evidence for establishment of dams in the Medieval as the factor causing divergence. Bayesian cluster analysis 
showed influence of stocked trout in several reference populations, but not in the focal lake and anadromous 
populations. Estimates of effective population size using a linkage disequilibrium method ranged from 244 to > 1,000 
in all but one anadromous population, but were lower (1 53 and 252) in the lake populations. 

Conclusions: We show that genetic divergence of lake-dwelling trout in two Danish lakes reflects establishment 
of water mills and impassable dams ca. 600-800 years ago rather than a natural genetic population structure. 
Although effective population sizes of the two lake populations are not critically low they may ultimately limit 
response to selection and thereby future adaptation. Our results demonstrate that populations may have been 
affected by anthropogenic disturbance over longer time scales than normally assumed. 

Keywords: Approximate Bayesian Computation, Divergence time, Effective population size, Habitat fragmentation, 
Isolation-with-gene-flow model, Microsatellite DNA 



Background 

Nearly all natural ecosystems, and the species that they 
encompass, experience anthropogenic pressure resulting 
particularly from habitat destruction, climate change, 
overharvesting and introduction of exogenous species 
[1-4]. Land use leading to habitat fragmentation is con- 
sidered one of the biggest threats to biodiversity [1], 
also encompassing short- and long-term negative gen- 
etic consequences [5-7]. Fragmentation problems are 
particularly severe in freshwater systems [8], where the 
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one-dimensional structure of river systems causes dams, 
weirs and other human-made obstructions to represent 
impassable barriers to many fishes and invertebrates. 
Accordingly, several studies have documented negative 
effects of dams on freshwater fish populations that de- 
pend on migration between different spawning, nursery 
and foraging habitats in different parts of their life 
cycle. In some cases this had led to extirpation of entire 
populations [9-12]. In other cases fragmentation due to 
dams has been shown to exert negative genetic effects 
caused by restricted gene flow and declining population 
sizes [13-17]. 

Populations within a species are often distributed across 
habitats that exhibit different environmental conditions. 
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These differences may lead to variation in local selection 
regimes and potentially local adaptation [18], the degree of 
which is determined by the opposing forces of local selec- 
tion on the one side and gene flow and random genetic 
drift on the other [19,20]. Whereas the limiting factors of 
gene flow and drift under pristine conditions are influ- 
enced by natural barriers in the landscape and habitat size 
and quality [21,22], anthropogenic habitat fragmentation 
leading to decreased gene flow and effective population 
sizes may substantially shift migration-selection-drift equi- 
libria [23]. The outcome for local adaptation is not trivial 
[24,25] . Reduced effective population sizes may cause drift 
to overwhelm directional selection [26] and limit evolu- 
tionary potential [27]. Conversely, reduced gene flow may 
at the same time increase possibilities for local adaption, 
as directional selection becomes a stronger evolutionary 
force than immigration [19,28], though in the longer 
term gene flow into small populations may benefit 
adaptive responses by introducing new variation. In 
total, when investigating adaptive divergence between 
populations in human-altered environments the ques- 
tion arises if local adaption reflects historical or more 
recent anthropogenically modified selection regimes 
and demographic parameters. 

Most salmonid fish species form local, genetically differ- 
entiated populations, many of which are locally adapted 
[29,30]. This is also the case for the brown trout {Salmo 
trutta), where several recent studies based on genetic 
markers and outlier scans [31-33], transcriptomics [34] 
and quantitative genetics experiments [35-38] have 
provided evidence for local adaptation. Specifically, in 
a common garden experiment a lake-dwelling brown 
trout population from Lake Hald, Central Jutland, 
Denmark (Figure 1) showed evidence for being adapted 
to higher water temperatures during incubation of eggs 
and larvae as compared to other populations [37]. In a 
second study, both the Lake Hald and another lake- 
dwelling trout population from Lake Moss0 (Figure 1) 
showed several outlier loci in microsatellite DNA out- 
lier scans encompassing anadromous, landlocked and 
hatchery trout populations [32]. 

Both Lake Hald and Lake Mosso are part of the same 
river system, the Gudena River (Figure 1). There are no 
natural barriers to migration within the system and 
historically anadromous trout would have access to all 
parts of the river. However, over the past ca. 800 years 
impassable dams have increasingly been established. 
During the Medieval water mills, often owned and man- 
aged by monasteries, were established at many rivers for 
manufacturing of textiles and grinding of flour. Whereas 
the first types of water mills built from ca. 1200 may not 
have completely and permanently blocked the rivers, 
from ca. 1400 new types of mills became established 
that involved permanent damming and created impassable 



barriers [39]. Lake Hald and Lake Mosso became isolated 
from the headwaters of the Gudena River system during 
that time and must be assumed to have been inaccessible 
to anadromous brown trout from the downstream part 
of the river system ever since. Hence, the question arises 
as to how historical anthropogenic fragmentation has 
affected the genetic structure, demographic parameters 
and potentially adaptive divergence of these popula- 
tions. Did the lake trout populations diverge naturally 
thousands of years ago, potentially as far back in time as 
the end of the last Glaciation ca. 12,000 years bp? Or is 
habitat fragmentation due to establishment of dams the 
major factor shaping the current genetic population 
structure? What are the genetic consequences of several 
centuries of reproductive isolation of the populations? 

In the present study we analyze variation at 44 micro- 
satellite loci in brown trout populations from Lake Hald, 
Lake Moss0 and anadromous trout from the lower part 
of the Gudena River (the Lillea River tributary). For refer- 
ence, we additionally include 10 other anadromous popu- 
lations and two hatchery strains that have been used for 
stocking in the region. Using individual Bayesian cluster- 
ing [40] we first assess if stocking has significantly affected 
the current Gudena River populations and genetic struc- 
ture. Next, using analyses based on Approximate Bayesian 
Computation (ABC) [41] and an Isolation-with-migration 
model [42] we estimate divergence time between anadro- 
mous trout from the lower Gudena River and the two lake 
populations. If the populations already were significantly 
reproductively isolated through natural processes thou- 
sands of years ago, we expect this to be reflected in the 
divergence time estimates. Conversely, if anthropogenic 
fragmentation is the primary determinant shaping the 
genetic structure of the populations, we expect estimates 
of divergence time to coincide with the establishment of 
dams ca. 500-800 years ago. Finally, based on results from 
the same analyses along with estimates of effective popula- 
tion size based on a linkage disequilibrium method [43] 
we assess the impact of dams on demographic parameters 
of the populations and discuss the consequences for 
extant and future adaptive responses. 

Methods 

Study populations 

Lake Hald (in the following abbreviated HAL); see 
Figure 1) covers an area of 3.3 km and is inhabited by 
brown trout that spawn in the tributaries. The outlet is 
Non Mollea River, the Danish name of which indicates 
that nuns from the nearby Asmild Monastery maintained 
a water mill at the river. Damming of the outlet and isola- 
tion of HAL from the Gudena River system is assumed to 
have taken place sometime in the period from 1400 to 
1500, although older mills have been constructed even 
earlier. Stocking of adult hatchery strain trout directly in 
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the lake took place in the 1970s-1980s, but a study based 
on allozymes suggested little if any genetic impact [44] . 

Lake Mosso (in the following abbreviated MOS; see 
Figure 1) covers an area of 16.9 km and is also inhab- 
ited by brown trout that spawn in the tributaries. The 
nearby 0m Monastery established a water mill and per- 
manent dam at Rye Ivtelle downstream of the lake 
around 1500. However, remnants have been found of 
older water mills constructed 100-200 years earlier. 
There has been limited stocking of hatchery strain 
trout in some of the tributaries, and transplants from 
the lower part of the Gudena River could potentially 
also have affected the genetic composition. 

The Lillea River (see Figure 1; abbreviated LIL) is a 
major tributary of the lower Gudena River. Access to the 
Gudena River for spawning anadromous trout was signifi- 
cantly reduced after construction of the Tange hydropower 
plant (see Figure 1) in 1920. The majority of the present 
anadromous trout population of the Gudena River now 
spawn in LIL downstream of the hydropower plant, and 
we considered trout from this river as representative of 
the trout population downstream HAL and MOS. The 
population could potentially be affected by stocked hatch- 
ery strain trout, in particular indirectly due to stocking 
elsewhere in the system. 

To investigate if stocked hatchery strain trout genetically 
affected the populations, which could influence estimation 
of divergence time and demographic parameters, we in- 
cluded two hatchery strains in the analyses, Vork Hatchery 
(VOR) and Harkser Hatchery (HAR). Stocking with hatch- 
ery strain trout has not been permitted in Denmark since 
2003, but was wide-spread previously. More than 80% of 
all stocked fish were derived from four strains that share a 
common history and have been shown to exhibit close 
genetic relationships [45]. VOR and HAR belong to this 
group of strains, which was originally founded by wild 
spawners from the Vejle and Kolding River on the East 
coast of Jutland. 

Finally, to assess possible sources of gene flow into the 
LIL population and to provide a comparison to estimates 
of demographic parameters and to populations known 
to be introgressed by stocked hatchery strain trout, we 
included samples from 6 anadromous trout populations 
on the Jutland West coast (Stora River, STO; Skjern 
River, SKJ; Varde River, VAR; Sneum River, SNE; Kongea 
River, KON; Ribe River, RIB), two populations from the 
Limfjord (Karup River, KAR; Skals River, SKA) and two 
populations from the Jutland East coast (Villestrup River, 
VIL; Kolding River, KOL) (see Figure 1). The West coast 
rivers, particularly SKJ and VAR, have previously been 
found to be strongly admixed with hatchery strain trout 
[46]. In contrast, KAR has been heavily stocked with 
hatchery trout, but nevertheless shows very low admixture 
[47]. SKA has not been stocked and to our knowledge the 



same applies to the VIL and KOL populations. However, 
KOL is known to be a main source for the founding of the 
quantitatively most important hatchery strains, including 
VOR and HAR, and could for that reason be expected to 
show close genetic relationships to these strains. 

All samples were collected by electrofishing conducted 
from 1999-2009 and encompassed multiple cohorts. Elec- 
trofishing was conducted by technical staff at the Danish 
Institute for Fisheries Research (now National Institute 
of Aquatic Resources, technical University of Denmark), 
who had all necessary permits from the Danish Ministry 
of Food, Agriculture and Fisheries for electrofishing and 
sampling of tissue and followed all required regulations. 
The fish were anaesthetized and a small piece of adipose 
fin was removed and stored in 96% ethanol, following 
which the fish were released. Sampled localities, sample 
abbreviations, ecotypes (lake-dwelling, anadromous, hatch- 
ery trout) and sample sizes are summarized in Table 1. 

Molecular markers 

The study was based on a subset of the data set from a 
previous publication [32]. This paper analyzed 74 micro- 
satellite markers and used outlier scans and landscape 
genomics approaches for detecting loci under possible 
selection. For the present study we omitted loci that 
were previously found to be outliers in terms of genetic 
differentiation, showed significant association with envir- 
onmental parameters, were known to be linked to func- 
tional loci, and/or showed distributions of allele sizes that 
would indicate strong deviation from a stepwise mutation 
pattern. In total, we analyzed 44 loci, which are listed in 
Additional file 1: Table SI. For technical details on DNA 
extraction and molecular analyses we refer to a previous 
publication [32]. 

Analysis of genetic variation and differentiation 

Exact tests for deviations from Hardy- Weinberg equi- 
librium, observed (H 0 ) and expected (H e ) heterozygos- 
ity and allelic richness was quantified with FSTAT 2.9.1 
[48]. Genetic differentiation between populations was 
quantified with an estimator of F ST [49] and significance 
tested by permuting alleles (10 4 times) among popula- 
tions using MSA 4.05 [50]. 

Bayesian cluster analysis 

Bayesian clustering implemented in STRUCTURE 2.3.4 
[51,52] was used for estimating the number of populations/ 
groups represented by the sampled individuals (k) and for 
estimating individual and population-level admixture pro- 
portions. The analysis assumed correlated allele frequen- 
cies and was furthermore based on the LOCPRIOR model 
[40], where information on the sample origin of individ- 
uals was used as prior information. For estimating the 
most likely k we conducted runs assuming k of 1 through 
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Table 1 Information on sampled populations 


Population 


Abbreviation 


Geographical region 


Ecotype 


Sample size 


Lake Hald 


HAL 


Eastern Jutland, Kattegat Sea (Gudena River system) 


Lake-dwelling 


32 


Lake Moss0 


MOS 


Eastern Jutland, Kattegat Sea (Gudena River system) 


Lake-dwelling 


32 


Lillea River 


LIL 


Eastern Jutland, Kattegat Sea (Gudena River system) 


Anadromous 


32 


Stora River 


STO 


Western Jutland, North Sea 


Anadromous 


34 


Skjern River 


SKJ 


Western Jutland, North Sea 


Anadromous 


53 


Varde River 


VAR 


Western Jutland, North Sea 


Anadromous 


36 


Sneum River 


SNE 


Western Jutland, North Sea 


Anadromous 


35 


Kongea River 


KON 


Western Jutland, North Sea 


Anadromous 


33 


Ribe River 


RIB 


Western Jutland, North Sea 


Anadromous 


31 


Skals River 


SKA 


Central Jutland, Limfjord 


Anadromous 


32 


Karup River 


KAR 


Central Jutland, Limfjord 


Anadromous 


32 


Villestrup River 


VI L 


Eastern Jutland, Kattegat Sea 


Anadromous 


32 


Kolding River 


KOL 


Eastern Jutland, Kattegat Sea 


Anadromous 


32 


Harkaer Hatchery 


HAR 


Hatchery 


Hatchery strain 


35 


Vork Hatchery 


VOR 


Hatchery 


Hatchery strain 


34 



16. Each run consisted of a burn-in of 10 5 MCMC steps, 
followed by 2x10 s steps. Ten replicates were conducted 
for each k. We plotted the probability of the data [(P(D)] 
and the ad hoc statistic AK (Evanno et al. 2005) which 
measures the steepest increase of the probability of k, 
using STRUCTURE HARVESTER [53]. CLUMPP [54] 
was used to find the consensus configuration based on the 
ten replicates for a given k, implementing the LargeK- 
Greedy algorithm. Finally, membership proportions of 
individuals and samples to the identified clusters were 
visualized using DISTRUCT [55]. 

Estimation of contemporary effective population size 

We estimated the contemporary effective population size 
of all populations using the linkage disequilibrium (or 
more precisely gametic phase disequilibrium) method 
[43,56]. For this purpose we used the software LDNE 
[57] and excluded alleles occurring at frequencies below 
0.05. Given that the samples analyzed included multiple 
age classes, we assume that the estimates are closer to 
representing N e (the effective population size per gener- 
ation) than Nb (the effective population size per cohort). 

We also tested for recent bottlenecks in MOS, HAL 
and LIL using the method and software BOTTLENECK 
[58,59] and assuming a two-phase mutation model with 
90% stepwise mutations [60]. 

Estimation of divergence time and demographic 
parameters 

We used Approximate Bayesian Computation (ABC) as 
implemented in DIYABC 2.0 [41] to assess posterior 
likelihoods of divergence time between LIL and HAL and 
between LIL and MOS, along with effective population 



sizes. It should be noted that DIYABC assumes that no 
gene flow occurs after populations have split. We assumed 
a generation time of 3.5 years [61]. We used wide and flat 
priors of N e for all populations [50; 10,000] and for diver- 
gence time, t (30; 3,000 generations, corresponding to 105 
and 10,500 years, respectively). We assumed a generalized 
stepwise mutation model [62] with a uniform prior distri- 
bution of mean mutation rate from 10~ 4 to 10 3 , a prior 
distribution of individual locus mutation rates from 10~ 5 
to 10~ following a Gamma distribution with mean deter- 
mined by the mean mutation rate across loci. The number 
of repeats per mutational event across loci was assumed 
to follow a geometric distribution with a uniform prior 
for the parameter P ranging from 0.1 to 0.6, whereas P 
for individual loci followed a Gamma distribution with 
mean determined by mean across loci and prior ranging 
from 10~ 2 to 9 x 10 1 . A number of summary statistics 
can be chosen for estimating posterior distributions of 
parameters, some of which are, however, partly redun- 
dant. We followed the approach of a previous study [63] 
and conducted three different sets of analyses for each 
scenario based on combinations of summary statistics 
found to be useful in previous studies; 1) mean number 
of alleles across loci within populations, mean expected 
heterozygosity within populations, mean value of M [64] 
within populations, F S t between populations and across 
loci, and index of individual assignment [65] across popu- 
lations [66]; 2) mean number of alleles across loci, mean 
expected heterozygosity, mean allele size variance within 
populations within populations, and the same statistics 
across populations [67]; 3) F ST between pairs of popula- 
tions, mean individual assignment likelihoods of individ- 
uals collected in one of a pair of populations and assigned 
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to the other [68], the mean number of alleles per locus be- 
tween pairs of populations, mean expected heterozygosity 
between pairs of populations and mean variance of allele 
size [69] between pairs of populations [70] . 

We pre-evaluated each scenario and set of priors by per- 
forming 10 4 simulations, conducting a Principal Compo- 
nent Analysis (PCA) based on the summary statistics and 
checking if the observed data set and the cloud of simu- 
lated data were congruent. The analyses were based on 
simulating 10 6 data sets, and the posterior distribution of 
parameters was estimated using the logit approach based 
on the 10 4 (1%) data sets closest to the observed data. Fol- 
lowing estimation of parameters we checked the model by 
simulating 10 3 data sets based on the predictive posterior 
distribution of parameters but using a new set of summary 
statistics, and conducting a PCA based on the summary 
statistics of these simulated data sets along with the 10 4 
data sets simulated based on the prior distribution (see 
above) and the summary statistics from the observed data. 
These model checks were conducted by switching sum- 
mary statistics among analyses. Hence, for analyses based 
on summary statistics set a) we checked the model using 
summary statistics set b), for the analyses using set b) we 
checked the model using set c), and for the analyses using 
set c) we checked the model using set a). 

We further estimated divergence time, effective popu- 
lation sizes and gene flow using the Bayesian, coales- 
cence and Markov Chain Monte Carlo (MCMC) based 
method IMa [42]. The method is based on an isolation- 
with-gene-flow model, where a single population back 
in time splits into two, which are subsequently con- 
nected by some gene flow. The following parameters 
are estimated: t, the point in time when the ancestral 
population split into two, qA, the effective population 
size in the ancestral population prior to splitting; qx and 
q 2 , the effective population size of population 1 and 2 after 
divergence; mi and rvi2, the migration rate from popula- 
tion 2 into population 1 and from population 1 into popu- 
lation 2. The parameters are scaled by mutation rate. To 
provide unsealed parameters we assumed the microsatel- 
lite mutation rate estimated using DIYABC (see above). 
Based on initial trial runs we found that a heated chain ap- 
proach encompassing 30 chains (parameters gl = 0.7 and 
g2 = 0.4), an initial burn-in of 2 x 10 6 MCMC steps 
followed by 2 x 10 6 steps yielded convergence. We as- 
sumed a generation time of 3.5 years (see above) and 
the following upper bounds on scaled priors: ql = q2 = 
10, qA = 100, ml = m2 = 175 and t = 2. We conducted 
three replicate runs for each population pair. 

Results 

Genetic variation and differentiation 

Summary statistics for all loci in all populations (ex- 
pected and observed heterozygosity, allelic richness, 



tests for Hardy- Weinberg equilibrium) are listed in 
Additional file 2: Table S2. Pairwise F ST ranged from 
0.002 (between the neighbouring populations RIB and 
KON) to 0.078 between HAL and the hatchery strain 
HAR. All pairwise F ST estimates were statistically sig- 
nificant except for that between RIB and KON. The two 
lake populations HAL and MOS were significantly dif- 
ferentiated from all other populations with F S t ranging 
from 0.025 to 0.062. Differentiation between HAL and 
the downstream anadromous LIL population within the 
Gudena River system was 0.033, whereas F S t between 
MOS and LIL was 0.037. Finally, F ST between HAL and 
MOS was 0.053 (Additional file 3: Table S3). 

Bayesian cluster analysis 

Analysis of the STRUCTURE output showed that AK 
was highest for k = 2, but with additional peaks for k = 4, 
7 and 12 (Additional file 4: Figure SI). At k= 2 one clus- 
ter was prevalent (>90%) in the two hatchery strains, 
VOR and HAR, whereas the other cluster was prevalent 
(>90%) in KAR, SKA, MOS and HAL. The other popula- 
tions showed significant proportions of both clusters 
(data not shown). Assuming k=7 provided the most 
detailed separation of populations (see below), whereas 
assuming /c=12 did not provide further biologically 
meaningful separation (data not shown). 

At k = 7, the hatchery strains VOR and HAR were 
characterized primarily by the "red" cluster in Figure 2, 
which also occurred at high frequency in the West coast 
populations, consistent with previous findings of ad- 
mixture with hatchery strain trout [46]. The cluster 
was virtually absent from KAR, SKA, VIL, HAL and 
MOS, but was present at low frequency (0.075) in LIL 
and high frequency (0.423) in KOL. This could reflect 
a common population history, as the two hatchery 
strains were founded by trout from the two East coast 
populations KOL and the neighbouring Vejle River, or 
it could represent admixture due to stocking. Compari- 
son of population and individual level cluster member- 
ship proportions (Figure 2) showed a distinct pattern 
for the stocked West coast populations, where some 
individuals showed strong admixture proportions of 
the "red" hatchery-specific cluster, whereas other were 
virtually non-admixed. In contrast, all individuals in 
LIL and KOL showed approximately equal proportions 
of the "red" cluster, indicating that a common population 
history rather than admixture underlies the results. We 
therefore assume that LIL is not admixed with hatchery 
trout and is representative of indigenous anadromous 
Gudena River trout. 

HAL and MOS were characterized by two different 
distinct clusters ("orange" and "yellow") whereas the two 
Limfjord populations KAR and SKA were characterized 
by a "green" cluster (Figure 2). The southern populations 
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Population Individual 
admixture admixture 



STO 

SKJ 

VAR 
SNE 
KON 
RIB 
KAR 
SKA 
VI L 
LIL 
KOL 
HAL 
MOS 
VOR 
HAR 



fc 




Anadr. W. coast 



Anadr. W. coast 



Anadr. W. coast 



Anadr. W. coast 



Anadr. W. coast 



Anadr. W. coast 



Anadr. Limfjord 



D Anadr. Limfjord 
Anadr. E. coast 




Anadr. E. coast 

Anadr. E. coast 

Lake 

Lake 

Hatchery 

Hatchery 



Figure 2 Admixture proportions estimated using STRUCTURE 
2.3.4 [51,52] and assuming k=7. The left panel shows population 
level admixture proportions whereas the right panel shows 
individual admixture proportions. 



on the Jutland West coast (SNE, KON, RIB) showed 
high proportions of a "dark blue" cluster that was near 
absent in other populations. Finally, a "light blue" cluster 
was found in high proportion in the East coast popula- 
tions VIL, LIL and KOL and to a smaller extent in some 



of the West coast populations. In total, the results suggest 
that the most likely sources of gene flow into LIL encom- 
pass other East coast populations and that the gene pool 
of this population, along with MOS and HAL are unlikely 
to be strongly affected by stocked hatchery strain trout. 

Contemporary effective population size 

Estimates of contemporary effective population size (N e ) 
using LDNE generally ranged from ca. 250 to > 1,000 in 
the anadromous trout populations, the exception being 
SKA, where the point estimate of Ne was 74 (Table 2). 
In RIB, N e was too high to be estimated (denoted by °°, 
meaning that sampling variance exceeded the signal 
from drift), but the lower 95% confidence interval was 
480. In four anadromous trout populations (including 
RIB) the upper 95% confidence interval could not be de- 
termined. Specifically for LIL, representing anadromous 
trout from the lower Gudena River system, N e was 288. 
A comparable N e estimate was obtained for MOS (252), 
whereas the estimate for HAL was lower (153). 

The BOTTLENECK tests [58] provided no evidence 
for recent bottlenecks in neither LIL, HAL nor MOS. 

Divergence time and demographic parameters 

The analyses using DIYABC yielded point estimates of 
mean mutation rate ranging from 2.90 x 10 4 to 3.85 x 
10 4 (Table 3). Estimates of divergence time for LIL-HAL 
ranged from 602 to 742 years, assuming different sets of 
summary statistics, whereas for LIL-MOS it ranged 
from 599 to 641 years, thus providing strong support 
for divergence caused by the establishment of water mills 

Table 2 Effective population size (N e ) estimates, obtained 
using the linkage disequilibrium method [43] 
implemented in the software LDNE [57] 



Population 


N e (95% CI) 


LIL 


288 (183-647) 


HAL 


153 (115-227) 


MOS 


252 (163-533) 


STO 


244 (166-445) 


SKJ 


259 (198-370) 


VAR 


245 (1 74-407) 


SNE 


429 (242-1706) 


KON 


1317 (358-oo) 


RIB 


oo (480-=) 


KAR 


537 (257-=) 


SKA 


74 (63-89) 


VIL 


369 (209-1401) 


KOL 


771 (307-°°) 


VOR 


91 (75-113) 


HAR 


248 (167-469) 
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Table 3 Summary of estimates of effective population size (W e ), divergence time and mutation rate (95% confidence 
intervals in parentheses) estimated using DIYABC [41], based on different sets of summary statistics as detailed in 
Methods 



Populations 


Summary 
statistics 


N e LIL 




N e HAL 


N e MOS 


Divergence time 
LIL-HAL (years) 


Divergence time 
LIL-MOS (years) 


Mutation rate 




L L-MOS 


Set 1 


7040 (316 


0-9630) 




3420 (922-8020) 




641 (258-1236) 


2.90x1 0" 4 
(1.63x10~ 4 -6.92x 


0" 4 ) 


L L-MOS 


Set 2 


5790 (216 


0-9430) 




3670 (1170-7420) 




620 (254-1386) 


3.85x1 0' 4 
(2.33x1 0~ 4 -8.12x 


0" 4 ) 


L L-MOS 


Set 3 


5930 (216 


0-9460) 




4070 (1180-8760) 




599 (217-1246) 


3.1 9x1 0" 4 
(1.76x10~ 4 -7.35x 


0" 4 ) 


L L-HAL 


Set 1 


7470 (3540-9730) 


1 980 (490-5820) 




602 (245-1 1 20) 




3.02x1 0" 4 
(1.68x10~ 4 -7.33x 


0~ 4 ) 


L L-HAL 


Set 2 


6830 (3010-9640) 


2900 (843-6510) 




742 (287-1715) 




3.71 X10" 4 
(2.25x1 0~ 4 -7.95x 


o~ 4 ) 


L L-HAL 


Set 3 


6800 (2920-9570) 


3470 (968-8360) 




686 (238-1425) 




3.18x10~ 4 
(1.75x10~ 4 -7.39x 


o~ 4 ) 



and associated impassable dams during the Medieval. Esti- 
mates of N e generated by DIYABC were substantially 
higher than those obtained using LDNE. Hence, point esti- 
mates of N e in LIL ranged from 5790 to 7470 across differ- 
ent analyses, in HAL it ranged from 1980-3470 and in 
MOS from 3420-4070 (Table 3). 

For the IMa analyses we assumed a mutation rate of 
3.0 x 10~ corresponding to the estimates provided by 
DIYABC. The results showed highly consistent outcomes 
across the three replicate runs for LIL-HAL and LIL- 
MOS, both concerning estimates of divergence time 
(Figure 3) and effective population size and migration 
rate (Table 4). The mean divergence time across the 
three runs was 789 years for LIL-HAL and 696 years for 
LIL-MOS (Table 4), corresponding well to the estimates 
obtained using DIYABC. There was virtually no statis- 
tical support for divergence having occurred naturally 
further back in time (> 1,500 years; Figure 3). Estimates 
of N e of the two isolated lake populations MOS and HAL 
were lower than for the downstream LIL population 
(Table 4). Average modes across three runs yielded esti- 
mates of 304 in MOS and 369 in HAL, whereas average 
mode of N e in LIL across six runs was 1585. Hence, the 
estimates were higher than those obtained using LDNE, 
but considerably lower than suggested by DIYABC. Esti- 
mates of N e of the ancestral population prior to divergence 
were high, generally on the order of 40,000 (Table 4). Gene 
flow estimates from LIL and upstream above dams to 
HAL and MOS were, as expected, virtually zero. Also, 
there was limited support for gene flow downstream from 
HAL to LIL, whereas low but non-negligible gene flow 
from MOS to LIL was indicated (Table 4). 

Discussion 

Divergence time estimates obtained using two different 
methods based on Approximate Bayesian Computation 



(DIYABC) and an Isolation-with-migration model (IMa) 
led us to reject the hypothesis that divergence of HAL 
and MOS from the anadromous LIL population is natural 
and occurred several thousand years ago. Instead, diver- 
gence time coincided remarkably well with establishment 
of water mills and impassable dams in the Medieval from 
ca. 1200 to 1500. There was less congruence regarding the 
demographic impact of the establishment of dams, as 
DIYABC consistently provided higher estimates of effect- 
ive population sizes than did IMa and estimates based on 
the linkage disequilibrium method (LDNE). 

Bayesian cluster analysis showed that the two lake 
populations HAL and MOS were genetically distinct, 
whereas the anadromous trout populations showed a 
regional genetic structure, reflecting ongoing gene flow 
and isolation-by-distance, as documented in other studies 
[32,46]. Importantly, these results also showed a minimal 
genetic contribution by stocked trout to the focal popu- 
lations LIL, MOS and HAL, which could otherwise 
complicate estimation of divergence time and demo- 
graphic parameters. 

In the following, we first discuss the reliability of the es- 
timates provided by DIYABC and IMa, and subsequently 
the impact of the impassable dams on genetic population 
structure and adaptive divergence. 

Reliability of estimates 

The conclusion that the two lake populations diverged 
from anadromous trout due to habitat fragmentation 
during the Medieval relies on the results using the 
methods IMa [42] and DIYABC [41]. In both cases we 
used wide and flat priors that should not strongly in- 
fluence posterior distributions. The mutation rate as- 
sumed, however, is a factor of uncertainty. For the IMa 
analyses we assumed a rate of 3.0 x 10" , estimated 
using DIYABC. This is quite close to estimates of 
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Figure 3 Probability density plots of divergence time estimates obtained using IMa [42]. A) Divergence time between LIL and HAL 
B) Divergence time between LIL and MOS. The results from three replicate runs are shown, indicated by different colours. 



mutation rate at dinucleotide loci in e.g. humans of 
2.73 x 1(T 4 [71] and 5.56 x 1CT 4 in common carp 
{Cyprinus carpio) [72]. Hence, the mutation rate as- 
sumed must be considered realistic, and it would require 



a much lower mutation rate (2.01 x 1(F 5 ) to change the 
estimate of divergence time for LIL-HAL from 789 years 
to, say, 12,000 years, coinciding with the end of the last 
Glaciation. 
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Table 4 Summary of results from three runs of IMa [42] with different starting points for each of the population pairs 
LIL-HAL and LIL-MOS 



Population 
pair 


Run 

number 


t 




<?r 






1a 


m, (LIL to HAL/MOS) 


m 2 (HAL/MOS to LIL) 


LIL-HAL 


a 


641 (339- 


-1366) 


1496 


(695-2403) 


497 (187-695) 


40625 (34125-50041) 


0.000 (0.000-0.008) 


0.003 (0.000-0.015) 




b 


804 (385- 


-1525) 


1588 


(895-2496) 


313 (196-678) 


41542 (33559-51041) 


0.000 (0.000-0.006) 


0.004 (0.002-0.013) 




c 


922 (431- 


-1528) 


1688 


(954-2602) 


296 (188-671) 


42292 (33958-51708) 


0.000 (0.000-0.006 


0.004 (0.002-0.014) 


LIL-MOS 


a 


688 (361- 


-1271) 


1588 


(1004-2313) 


321 (171-504) 


39376 (32874-48042) 


0.000 (0.000-0.002) 


0.011 (0.006-0.020) 




b 


712 (269- 


-1317) 


1529 


(895-2346) 


287 (146-537) 


39874 (32042-49208) 


0.000 (0.000-0.002) 


0.012 (0.004-0.021) 




c 


687 (360- 


-1225) 


1622 (1012-2287) 


304 (170-506) 


40625 (33542-48707) 


0.000 (0.000-0.002) 


0.011 (0.006-0.020) 



f Denotes divergence time {in years) between the two populations assuming a standard mutation rate of 5.56x1 whereas t est mut denotes divergence time 
(in years) assuming a mutation rate of 3.00x1(T 4 , reflecting the mutation rate estimates obtained using DIYABC. q, denotes effective population size of LIL, q 2 
effective population size of HAL or MOS and q A effective population size of the ancestral source population prior to divergence, m, denotes migration rate per 
generation from MOS or HAL into LIL, and m 2 denotes migration rate from LIL to HAL or MOS. The analyses are based on coalescence, but to increase clarity 
estimates of m, and m 2 are reported in forward-time. All estimates consist of modes and 90% credible intervals {in parentheses). 



Whereas the congruence of results from different rep- 
licated runs of IMa suggests good convergence, a poten- 
tial problem consists in the simplified model assumed, 
where LIL exchanges migrants with either HAL or 
MOS. In reality, both populations could simultaneously 
contribute to gene flow, but even more importantly LIL 
would be subject to gene flow from other anadromous 
populations. We are unable to resolve the magnitude of 
this potential problem, but we note that the results ob- 
tained generally make sense. For both lake populations 
the estimate of gene flow from LIL was virtually zero, as 
would be expected given the problems of passing the 
dams upstream. In contrast, there was a signal of down- 
stream gene flow into LIL from MOS. Gene flow in this 
direction is certainly a realistic possibility, as fish can 
pass the dams downstream by simply being flushed with 
the current. 

N e estimates of HAL and MOS obtained using IMa 
were of the same order of magnitude, albeit higher than 
estimates obtained independently using LDNE. In con- 
trast, estimates of N e in LIL obtained by IMa were more 
than 5 times times higher than those derived from 
LDNE. Although recent bottlenecks could explain this 
result, there was no evidence for this to have taken 
place. Alternatively, the result could reflect gene flow 
from other anadromous populations into LIL, thereby 
inflating historical N e estimates (as in IMa) compared to 
contemporary estimates (as in LDNE) that pertain to 
the current or only a few generations back in time. The 
complexities of estimating N e in populations showing 
some degree of geographical continuity are increasingly 
acknowledged [73], and we tentatively suggest that the 
N e estimate in LIL based on IMa should be interpreted 
more broadly as encompassing LIL and neighboring 
populations contributing to gene flow. 

The DIYABC analyses appeared robust towards the 
choice of summary statistics, as evidenced by similar 
outcomes of the different analyses. However, unlike IMa 



it is assumed that there is no gene flow following splitting 
of populations. In the present case, this is probably not a 
major problem, as IMa suggested low gene flow between 
LIL and HAL/MOS, and divergence time estimates 
were quite similar between methods. Effective popula- 
tion size estimates, however, differed strongly between 
IMa and DIYABC, with the latter providing 5-10 times 
higher estimates than IMa. In both cases the estimates 
are "historical", i.e. encompassing many generations, al- 
though IMa estimates N e from the ancestral population 
prior to divergence and N e of the separate populations 
after divergence, the latter thereby giving more weight 
to the recent past (in the present case the centuries 
since establishment of dams). Most summary statistics 
applied in DIYABC, such as expected heterozygosity, 
Fst an d numbers of alleles should primarily reflect his- 
torical N e over a very long time scale [74], with the ex- 
ception of M [64] which detects recent bottlenecks. 
Differences of the time scale over which N e is estimated 
could therefore explain the discrepancies of results. We 
emphasize the N e estimates obtained by IMa as being 
most realistic and relevant in a conservation context; 
they are generally congruent with estimates of contem- 
porary N e obtained using LDNE and with estimates 
from Danish brown trout populations obtained using 
temporal methods, that for anadromous populations in 
Danish rivers generally range from ca. 250 to > 1,000 
for large river systems [46,61]. 

In total, we find the estimates of divergence time ob- 
tained by IMa and DIYABC to be robust, whereas N e 
estimates differ, probably reflecting the different time 
scales that they apply to. 

Impact on genetic population structure and adaptive 
divergence 

Our results suggest that the two lake populations HAL 
and MOS were previously part of a coherent anadromous 
trout population inhabiting the Gudena River system. The 
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Gudena River is the largest watershed in Denmark, and 
the effective population size of the ancestral trout 
population prior to fragmentation has presumably been 
high; > 1,000, given the estimates obtained from the 
KON and RIB populations (see Table 2 and [46]). Gene 
flow has occurred to and from other such populations 
enabling introduction of new adaptive variation and 
potentially both enhancing and limiting local adaptation, 
depending on local selection regimes and migration-drift- 
selection equilibrium [24,25]. Following the establishment 
of water mills, HAL and MOS have been disconnected 
from these dynamics for several centuries. 

What is the primary impact of this fragmentation? 
Contemporary N e was estimated to 153 and 252 for HAL 
and MOS, respectively. Although these values cannot be 
considered high, they should also not be cause for imme- 
diate conservation concern. Indeed, they are higher than 
several estimates of N e in undisturbed resident brown 
trout populations from Sweden [75,76] that have been iso- 
lated naturally presumably for even longer time spans that 
the MOS and HAL populations, and are in fact compar- 
able to N e in some of the anadromous populations in- 
cluded in the study (Table 2). It should be considered, 
however, that several of these latter populations included 
for reference have likely declined recently, as documented 
by comparison of historical and contemporary samples 
from the populations [46] . 

We suggest that the most important consequence of 
the dams concerns local adaptation and evolutionary 
potential. Even if trout ascend the dams downstream 
and migrate to sea, they will be unable to return to 
spawn. This could impose strong selection against ana- 
dromy, although the response to selection depends on 
heritability of the trait. There are few studies available 
that have estimated heritability of anadromy in salmonids. 
However, a study based on pedigreeing a wild brook trout 
(Salvelinus fontinalis) population found heritability of life- 
history tactics (anadromy versus residency) as high as 
0.52-0.56 [77]. If this is also the case in brown trout, then 
the 600-700 years since establishment of dams corre- 
sponding to 170-200 generations should have left ample 
opportunities for selective responses to occur. On the 
other side, a recent study of trout in HAL documented 
that 15% of individuals aged between 1 and 3 years and 
with a length exceeding 12 cm left the lake and would 
potentially undertake migration to the sea, whereas 40% 
migrated into the lake and 44% remained in the tribu- 
taries [78]. Hence, a potential for long-distance migra- 
tion involving anadromy seems still to be present in the 
population, although it is unknown if a larger propor- 
tion of individuals would have left the lake prior to the 
establishment of dams. 

The lack of immigration from other populations into 
the lakes should in the short term reduce influx of 



locally maladaptive alleles thereby shifting the migration- 
selection balance in favor of local adaptation. A common 
garden experiment including both LIL and HAL trout 
demonstrated significantly different temperature-related 
reaction norms for early life history traits, with HAL 
showing adaptation to higher incubation temperatures 
during winter owing to the spawning tributaries being 
fed by ground water [37]. Whereas this selection regime 
would have been similar prior to the establishment of 
dams, local adaptation would be expected to be reduced 
depending on the rate of immigration from other popu- 
lations and their degree of maladaption [79] . 

In the long term, it is expected that reproductive isola- 
tion of the lakes combined with the relatively low effective 
population sizes should impose limits to the response to 
selection, due to the fact that loss of potentially adaptive 
variation is not counteracted by introduction of new vari- 
ation through migration. This is all the more serious, as 
anthropogenic pressure increases the need for populations 
to adapt to environmental change [2,4,80], and specifically 
for brown trout it is predicted that its future distribution 
will become significandy reduced due to loss of suitable 
habitat [81]. The requirements for maintaining evolu- 
tionary potential have been expressed in the classical 
"500 rule" [27] (but see also [82] and [83] for discus- 
sion), stating that N e should be 500 or more for new 
mutations to balance loss of quantitative variation by 
drift. This criterion is not fulfilled in HAL and MOS. 
Assuming no mutation, the limit to response to selec- 
tion, R(°o) should be given by the expression R(°°) = 2N e 
R(l), where R(l) denotes the initial response to selec- 
tion. Moreover, the time until 50% of the response, tso% 
is roughly equal to lAN e [84]. If we assume that HAL 
became isolated 700 years ago, corresponding to 200 
generations and that N e is 153, then tso% is 214 genera- 
tions, approximately equal to the number of generations 
that have elapsed since isolation of the population. Al- 
though this is based on simplified assumptions, it never- 
theless suggests that there should still be adaptive 
potential within HAL and the presumably larger MOS 
population. Of course, different conclusions could apply 
to other anthropogenically isolated populations with 
much lower N e . 

Conclusions 

Landscapes are presently being strongly modified by 
humans, but in some regions, such as large parts of 
Europe this is a process that has already been ongoing 
for centuries and even millennia. Hence, it is often difficult 
to ascertain that a given genetic population structure is 
natural rather than the result of anthropogenic modifica- 
tions. In this study, we show that genetic divergence of 
lake-dwelling trout in two Danish lakes reflects establish- 
ment of water mills and impassable dams ca. 600 years 
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ago rather than a natural genetic population structure. 
Also, the populations have historically been part of a larger 
system of anadromous brown trout. The reproductive iso- 
lation of the lakes is likely to have affected adaptive diver- 
gence among populations, in the short term by shifting 
the migration-selection balance. However, in the long term 
lack of gene flow combined with relatively low effective 
population size may compromise evolutionary potential. 

The results demonstrate that caution is warranted 
when analyzing genetic population structure and inter- 
preting it as "natural". Even though anthropogenic modi- 
fication of habitats has accelerated during the past 
century, it may nevertheless be a process that has been 
ongoing over a much longer time scale with significant 
consequences for the affected populations. Our results 
exemplify the use of molecular dating for assessing an- 
thropogenic factors affecting wild populations. This is 
similar to other recent studies demonstrating that popu- 
lation decline in Borneo orang-utans (Pongo pygmaeus) 
is more recent and drastic than previously assumed [85], 
and conversely that Californian fishers (Martes pennanti) 
became fragmented and experienced population declines 
prior to European setdement [86]. 
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